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An ubiquitous mechanism for water-like anomalies 
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PACS 61 . 20 . Ja - Computer simulation of liquid structure 

PACS 61 . 25 . Em - Molecular liquids 

PACS 65 . 20 . -w - Thermal properties of liquids 

Abstract. - Using collision driven molecular dynamics a system of spherical particles interacting 
through an effective two length scales potential is studied. The potential can be tuned by means of 
a single parameter, A, from a ramp (A = 0.5) to a square-shoulder potential (A = 1.0) representing a 
family of two length scales potential in which the shortest interaction distance has higher potential 
energy than the largest interaction distance. For all the potentials, ranging between the ramp and 
the square-shoulder, density and structural anomalies were found, while the diffusion anomaly 
is found in all but in the square-shoulder potential. The presence anomalies in square-shoulder 
potential, not observed in previous simulations, confirm the assumption that the two length scales 
potential is an ubiquitous ingredient for a system to exhibit water-like anomalies. 
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' Introduction. — Some liquids are Icnown as anoma- 
lous liquids since they exhibit unexpected behavior upon 
variations of its thermodynamic conditions. Water is the 
canonical example of those anomalous liquids. 
I Water expands upon cooling at fixed pressure [1], dif- 
fuses faster upon compression at fixed temperature [2,3], 
and becomes less organized upon increasing density - or 
equivalently upon compression - at constant tempera- 
ture [4]. These are the density, diffusion, and structural 
anomalies. 

I The region where these anomalies occur form nested 
domes in the density-temperature diagram [4] - or 
pressure-temperature diagram [5]. The structural 
anomaly domain occupies the outer region of the pressure- 
temperature phase-diagram and the density anomaly re- 
gion is the innermost region. The diffusion anomaly region 
lies between these two domains [4,5]. This is the hierarchy 
of anomalies of water. 

Water-like anomalies are also found in other liquids. For 
example, density anomaly was found experimentally in liq- 
uid Te [6], S [7,8], and GeisTegs [9]. Simulations for silica 
[10], silicon [11], and liquid beryllium [12] show that den- 
sity anomaly is also present in these materials. Diffusion 
and structural anomalies were found for silica [10, 13-15] 
and silicon [16] and structural anomaly is reported for liq- 
uid beryllium [12, 17] through simulations. 

In water, the density anomaly is due to the hydro- 



gen bonds. The compression of a hydrogen-bonded local 
environment leads to an increase in entropy, or, equiv- 
alently, that a local hydrogen-bonded environment pos- 
sesses a lower density than a non-bonded system would 
exhibit [18]. However, there is no hydrogen bond in Te, S 
and GeisTeg.g. 

Therefore, how can the anomalies be explained for bond- 
ing and non-bonding systems? In order to address this 
question, instead of looking for the specific mechanism be- 
hind the density anomaly in water or in silica, one has to 
find the universality behind it. 

The simplest framework in which one can look into the 
physics of anomalies is given by the two length scales po- 
tentials. The idea is that in principle an anisotropic in- 
terparticle potential can be modeled as effective poten- 
tial [19]. We will examine one class of effective two length 
scales potential in which a compression-based competition 
arises between particles population in the second and first 
shells. 

This assumption was confirmed in a number of isotropic 
two scales potentials [20-31] in which density, diffusion, 
and structural anomalies were found. This was also shown 
to be correct in anisotropic potentials in which the two 
length scales emerge from the mapping of the anisotropic 
potential in a equivalent spherical symmetric potential 
[19,32]. Ramp-like potentials have demonstrate to be par- 
ticularly useful since they can describe the effective in- 
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Fig. 1: Interparticle potential studied in this work. The poten- 
tial can be tuned by means of the parameter A, ranging from 
a ramp (A = 0.5) to a square-shoulder potential (A = 1.0). See 
the text for more details. 




teraction between clusters of water molecules [33,34]. In 
fact, Yan at al. showed a quantitative agreement between 
the phase-diagram for the ramp potential (A = 0.5 in Fig. 
1) and that one for the TIP5P molecular model for wa- 
ter [34]. The authors show that a central water molecule 
interacting with its four nearest neighbors can be modeled 
effectively through a ramp potential. 

A ramp liquid was also used to mimic water in a system 
composed by a mixture of water and hard-sphere particles 
[35] . The solvation thermodynamics of the ramp and hard- 
sphere mixture describes well the qualitative behavior of 
the water-like solvation thermodynamics [35]. 

For studying the water-like anomalies the exceptional 
case is the square-shoulder potential in which no anomalies 
were reported yet (Fig. 1 with A — 1.0) [36]. This led to 
the idea that anomalies are present in ramp- like potentials 
(Fig. 1 with A = 0.5) but not in shoulder-like potentials. 

The aim of this paper is to propose that two scales po- 
tentials, potentials in which two preferred distances are 
present, exhibit water-like anomalies. It will be shown 
that in some cases the anomalies are in an inaccessible 
region, as inside a crystal phase. This is the case for the 
square-shoulder potential [36] . 

The Model. — In order to test our assumption we de- 
veloped a tunable potential ranging from a ramp potential 
to a square-shoulder potential. The potential is given by 



Uir) 



oo, r < ao 

h{r), ao <r < ai 

i'2{r), CTi < r < 0-2 
0, (72 < r. 



(1) 



where 01 (r) [Uo {<7i -r) -Ui (do - r)] / (cti - do) and 
(j)2{r) = Ui (o-2 ~r) / ((72 " o-i). 

A single parameter A is used to tune the potential from 
a ramp (A = 0.5) to a square-shoulder (A = 1.0) where 
ci = Co + A ((72 — ^o) and Ui = XUo- 

In this work 02/ Uo = 1-75 and the potential was approx- 
imated by discrete steps in the same spirit of Ref . [37] , in 
such a way that the discrete energy-step is AJ7 = 0.025?7o. 



Fig. 2: Pressure-temperature (P-T) diagram for each A-case. 
Each line in the P-T diagram corresponds to an isochore: in 
the upper panel (A = 0.5), from bottom to top, they are p* = 
0.20, . . . , 0.34, and 0.35. In the middle (A = 0.6) isochores 
are p* = 0.19, 0.20, 0.21, 0.22, 0.23, 0.232, 0.24, . . . , 0.33, and 
0.34. Finally, in the lower panel (A = 1.0) p* = 0.17, 
0.29, and 0.30. The solid bold line connects the temperature 
of maximum density (TMD) points and it shrinks to a small 
region for A = 1.0. See the text for details. 



Systems with 0.5 < A < 1.0 were analyzed for density, dif- 
fusion, and structural anomalies. 

Using the collision driven molecular dynamics tech- 
niques [38] , systems with N = 500 identical spherical par- 
ticles of mass TO, interacting through the potential Eq. 
(1) and A — 0.5,0.6,1.0 were studied. These particles 
were confined into a cubic box with volume V and pe- 
riodic boundary conditions. The equilibration and pro- 
duction times were 500 and 1000 respectively, in units of 
Co^/ m/Uo (time units). The rescaling of the velocities 
scheme was used for every 2 time units in order to reach 
and keep the desired temperature. 

Simulations for locating the anomalies at the pressure- 
temperature phase-diagram were initialized with the sys- 
tem in the fluid phase. This procedure makes possible to 
sample the metastable liquid phase inside the solid phase. 

For estimating the melting line the system was initial- 
ized with particles in a face centered cubic configuration. 
After 500 time units it was checked if the system remains 
solid or if it has melted. This was done by checking the 
diffusion coefficient and the shape of the pair distribution 
function. 

This process gave an estimate of the melting line in the 
pressure-temperature phase-diagram in excellent agree- 
ment with the one presented in Ref. [37] for the case in 
which A = 0.5. 

Pressure, P, was calculated by means of virial [39] and 
diffusion, D, was derived from the mean square displace- 
ment [39]. The translational order parameter, t, was cal- 
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Fig. 3: Difltusion coefficient against density with each line cor- 
responding to an isotherm. Points are simulated data and 
continuous lines are fifth order polynomial fit from the data. 
Dashed lines connect maxima and minima in the D(p, T — 
constant) functions and they have the same meaning as in Fig. 
6. Between these extrema, diffusion anomalously increases un- 
der increasing density - just like water does [2] . In the upper 
panel (A = 0.5 case) the isotherms (from bottom to top) are 
T* = 0.08, 0.09, 0.14, and 0.15. In the middle panel 

(A = 0.6 case) the isotherms are T* = 0.08, 0.09, . . . , 0.13, and 
0.14. Finally, in the lower panel (A — 1.0 case), the isotherms 
are T* = 0.11, 0.12, 0.13, 0.14, and 0.15. 




Fig. 4; Translational order parameter, t, versus density for 
fixed temperatures. Points are simulated data and lines con- 
necting the points are fifth order polynomial fit from the data. 
Dotted lines bound the region where t decreases under increas- 
ing density and they have the same meaning as in Fig. 6. 
In the upper panel (A = 0.5) we show (from top to bottom) 
the temperatures T* = 0.08, 0.09, . . . , 0.14, 0.15, 0.17, 0.19, 
. . . , 0.33, and 0.35. In the middle panel (A = 0.6) are shown 
r* = 0.08, 0.09, . . . , 0.21, and 0.22. Finally, in the lower panel 
(A = 1.0) temperatures T* = 0.11, 0.12, 0.15, 0.16, 0.18, 
0.20, 0.22, and 0.24 are shown. 



Fig. 5: Orientational order parameter, Qg, versus density 
for fixed temperatures. Points are simulated data and lines 
connecting the points are fifth order polynomial fit from the 
data. Dashed-dotted lines connect the maximum values for 
Qq{p,T = constant) and they have the same meaning as in 
Fig. 6. The temperatures shown in the three panels are the 
same as in the Fig. 4. 



culatcd as [4] 



t = 



1.9(0 - 



(2) 



where ^ = rp^^^ is the interparticle distance r divided by 
the mean separation between pairs of particles p~^l^ . g{(^ 
is the pair distribution function and is a cut-off distance. 
In this work was used t,c = where L = V^^"^. For 

a completely uncorrelated system (ideal gas) g ~ 1 and t 
vanishes. In a crystal, a translational long-order {g ^ 1) 
persists over long distances making t large. 

The orientational order parameter [40], Qq, was com- 
puted using the strategy introduced by Yan el. al [29]. 
Qe was calculated as follows. First the spherical harmon- 
ics, YJ^, associated to each particle i and its k neighbors, 
were obtained by 



Ira 1 Yij > 



(3) 



The 



' stands for the average over the k neighbors. 
Details on the calculation of can be found in [21,29]. 
The orientational order parameter associated to each par- 
ticle i is then given by summing the spherical harmonic 
over all orders m for a fixed degree £ namely 



47r 



2£-t-l ^ 



Iml 



1/2 



(4) 



Then it was chosen i ~ & for characterizing the local 
order [21,29,32], so the orientational order parameter be- 
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comes 



1 ^ 



(5) 



that is the mean value of Qq over all particles of the sys- 
tem. 

The parameter Qg assumes its maximum value for a 
perfect crystal and decreases as the system becomes less 
structured. For a completely uncorrelated system (ideal 
gas) Q^Q = 1/y/k. In this work fe = 12 neighbors. For a 
crystal, the Qe value depends on the specific crystalline 
arrangement and on the number of neighbors taken into 
account. 

Pressure, diffusion, temperature, and density, p ~ N/V, 
are given in reduced units as P* = Pa^/Uo, D* = 



D {m/Uof' loo, T* = fcsT/[/o, and p* 



The Results. - The potentials with 0.5 < A < 1.0 
were tested for the presence of: a minimum at the iso- 
chores in the pressure-temperature phase-diagram; a max- 
imum and a minimum in the diffusion constant as a func- 
tion of density at constant temperature diagram, a maxi- 
mum and a minimum in the translational order parameter 
versus density at constant temperature, and a maximum 
in the orientational order parameter against density at 
constant temperature. Then the position of the anomalies 
were compared with the melting line. 

Fig. 2 illustrates the pressure-temperature (P-T) phase- 
diagram for the potential Eq. (1) with A — 0.5, 0.6 and 
1.0. Each line corresponds to an isochore (see the figure 
caption for details). Some isochores have minima, which 
define the temperature of maximum density (TMD). For 
the pressures and the temperatures inside the TMD, the 
system expands upon cooling under fixed pressure, thus 
this region is known as the density anomaly region [37]. 
As A increases, the density anomaly region shrinks and 
reduces to a very small region that we identify with a 
simple point the A = 1.0 case. This reduction of the den- 
sity anomalous region is consistent with simulations for 
0.6 < A < 1.0 (not shown) 

Fig. 3 shows the diffusion constant versus density at 
fixed temperature for A = 0.5, 0.6 and 1.0. For the A < 0.7 
cases (not shown) the diffusion constant versus density for 
a fixed temperature has a local maximum and a local min- 
imum (represented by a dashed line in Fig. 3). Between 
these local extrema the diffusion coefficient increases upon 
increasing density and this region is known as the diffusion 
anomaly region. 

Figs. 4 and 5 illustrate the translational and the orien- 
tational order parameters versus density for constant tem- 
peratures for A = 0.5, 0.6, and 1.0. The translational order 
parameter has a local maximum and a local minimum for 
an interval of temperatures. Between these local extrema, 
the parameter t decreases under increasing density. An 
anomalous t parameter was observed for all A-cases in the 
stable liquid phase, even for A = 1.0 where no anomalies 
were reported before. The t anomalous region is bounded 
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Fig. 6: Interparticle potential studied in this work with A = 
0.5,0.6 and 1.0 (left) and tlie corresponding results we have 
found (right). Dotted lines enclose the region where t decreases 
upon increasing density and dashed-dotted line mark the max- 
imum in the Qe parameter. The region between the dashed- 
dotted line and the high pressure dotted line is known as the 
structural anomaly region (see the text for details). Dashed 
lines bound the region of diffusion anomaly, and solid line de- 
termine the region where density anomaly occurs. For the case 
with A = 1.0 the density anomaly line shrinks into a very small 
region and the diffusion anomaly region lies inside the solid 
phase, becoming inaccessible. The shadowed region is bounded 
by the estimate of the melting line. 



by the dotted lines in Fig. 4. 

The orientational order parameter, Qe, has a maximum 
in the Qg ^ P plane for each isotherm as shown in Fig. 
5. For densities higher than the density of maximum Qg 
the orientational order parameter decreases as the density 
increases - the same trend seen for t. 

For each temperature the density of maximum Qe fies 
between between the density of minimum t and the den- 
sity of maximum t. Consequently for densities between 
the density of maximum Qg and the minimum t is the 
structural anomaly region. This region is illustrated in 
Fig. 6. 

Fig. 6 summarizes our findings for A = 0.5, 0.6, and 1.0. 
It shows the regions in the P-T phase-diagram where the 
density, diffusion, and structural anomalies are located. It 
also show the location of the melting line that bounds the 
region where the system becomes solid. 

For A = 0.5 - the ramp potential case - the high pressure 
dotted line and the dashed-dotted line bound the region 
of structural anomalies, where the translational and orien- 
tarional order parameters decrease with density. A dashed 
line encloses the diffusion anomaly region where particles 
diffuse faster upon increasing density. The solid line con- 
nects the temperatures of maximum density limiting the 
density anomaly region. The border of the shadowed re- 
gion estimates the outer limit for the melting line. For 
the ramp potential the region of density, diffusion, and 
structural anomalies are in the stable region of the pres- 



sure temperature phase-diagram. All these anomalies for 
the ramp case are well documented [32,37] and a further 
discussion on these results is unnecessary. 

For the potentials with A = 0.6 and 1.0 shown in Fig. 
6 the lines and shadowed region have the same meaning 
as for the case with A ~ 0.5. Comparing our results for 
A = 0.5, 0.6, and 1.0 an interesting trend is observed. 
Once the potential become "harder" - going from a ramp 
to a square-shoulder - three effects arc quite evident. 

The first effect of the change of A is related to the mo- 
bility of particles in the system. The diffusion anomaly 
region shrinks as A increases, moving to lower tempera- 
tures becoming inaccessible for the case in which A = 1.0. 
This is in accordance with the results obtained by Netz 
et al. [41]. These authors have observed that as the dis- 
cretization of the ramp potential becomes less coarser, 
with corresponding increasing in the energy steps, the 
diffusion anomaly region shrinks and migrates to lower 
temperatures into the density anomaly region. Lattice 
models exhibit this same effect since the lattice structure 
plays the role of a coarsely discretized energy barriers am- 
bient [42,43]. In both cases the lines in the pressure tem- 
perature phase-diagram defining the border between the 
density and diffusion anomalous regions cross for a certain 
choice of parameters what is also observed for A = 0.6. 

The second effect is related to the melting line. As A 
increases the estimate of the melting line goes towards high 
temperatures, engulfing the region of water- like anomalies. 
This could explain why no anomalies were reported for the 
case with A = 1.0. 

The third effect is related to both the structural and 
density anomaly regions. As A increases the temperature 
of these regions are not drastically affected. This is con- 
sistent with the results obtained by Netz et al. [41] , where 
the authors have observed that the discretization of the 
ramp potential does not affect the thermodynamics. In 
the current A increases the potential not only be- 

comes more discretized but also there is a change in the 
potential energy associated with each scale. This leads to 
a shrink in the pressure range of the anomalous region. 

We can also gain some insight on the relationship be- 
tween structure of the system, shape of the effective inter- 
particle potential, and presence of anomalies by analyzing 
the order map, i.e., the t — Qq plane [4,21,29,32]. Fig. 7 
show our results. 

The paths formed by the points in the order map for wa- 
ter collapse into a single line inside the structural anomaly 
region. This means that in the water case the order param- 
eters are strongly coupled. For silica [10] and beryllium 
fluoride [17], also anomalous tetrahedral liquids, the ori- 
entational and translational order parameters are weakly 
coupled, since they develop a two-dimensional region in 
the order map of such liquids. For all A considered in this 
work we have a silica- and beryllium fluoride-like behav- 
ior for the paths in the order map or our model. This 
means that the two scales potentials are hybrid models, 
in the sense that they can exhibit water-like hierarchy of 
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Fig. 7: Order map for the potential Eq. (1) with A = 0.5, 0.6, 
and 1.0. The arrows indicate the direction of increasing density 
for fixed temperatures. In the upper panel, the temperatures 
are (from top to bottom) T* = 0.08, . . . , 0.14, 0.15, 0.27, 0.33, 
and 0.35. In the middle panel, T* = 0.08, . . . , 0.13, 0.14, 0.18, 
and 0.22. Finally, in the lower panel, T* = 0.11, 0.14, 
0.15, 0.18, 0.20, and 0.22. 

anomalies [21,32] but reproduce the order map of silica 
and beryllium fluoride instead of water. The understand- 
ing of the coupling mechanism of the order parameters is 
not clear in the literature, and wc believe this could shade 
some light into important aspects of anomalous fluids. 

Next, a remarkable feature of the order map shown 
in Fig. 7 is that the inaccessible region is virtually A- 
independent. In all panels, the inaccessible region is 
bounded by a straight line given hy t ~ axQ^ + b\, with 
ao.5 = 12.05, ao.6 = 12.20, and ai.o = 12.92; 60.5 = -2.70, 
60.6 = —2.72, and 61.0 = —2.84. The difference between 
the extreme values of a\ and b\ quantities is less than 
7.5%. This means that the order in all A-systems respond 
roughly constant upon compression. In this sense we are 
lead to believe that the two scales feature, which is present 
in all A-cases considered here, is most important than en- 
ergetic barriers differences between them towards the ap- 
pearance of anomalies. Indeed, it was shown that struc- 
ture and water-like anomalies can be linked by means of 
the excess entropy [44] , stressing the importance of the role 
played by structural parameters into the understanding of 
the water-like anomalies. 

In resume, the three anomalous regions respond differ- 
ently to the change in the potential. The diffusion anoma- 
lous region shrinks in the pressure temperature phase- 
diagram and disappears in T* < 0.08 for A = 0.65 (not 
shown); the density anomalous region shrinks in pressure 
and reduces to a small region in A = 1.0; the structural 
anomalous region also shrinks in pressure but it is still 
in the stable region of the pressure-temperature phase- 
diagram for A = 1.0. The order map analysis show that 
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the inaccessible region in the t — Qq plane is virtually in- 
dependent of A and the order parameters are uncoupled, 
differently of water but similar to silica and beryllium flu- 
oride. 

Conclusions. — In this paper we have shown that 
two scales effective potentials always reproduce water-like 
anomalies. In this sense, any liquid material in which this 
kind of effective interaction is present is able to be an 
anomalous liquid. In some cases the anomalous regions are 
located in the pressure-temperature phase-diagram inside 
the region where the solid phase is the most stable, being 
inaccessible either with experiments or with equilibrium 
simulations. 

In water the two scales distances ~ two energetic- 
competing preferable distances - clearly arise from the for- 
mation and breaking of hydrogen bonds. In other anoma- 
lous liquids this process come from different competing 
interactions but the effective final interparticle potential 
must be a two scales potential. 

We believe that the knowledge of this mechanism can 
be of great interest for industries. The domain of this 
mechanism could lead to the development of new materi- 
als - polymers, for example - in which some anomalous 
properties could be used in the manufacture of substances 
close to its supercooled region. A polymer in which its 
molecules diffuses faster upon compression is an example 
of a possible application of these findings. 
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